#!/usr/bin/env Rscript

# Unit tests for 3D multitype point pattern functionality

library(testthat)
library(spatstat.geom)
library(spatstat.explore)

context("3D Multitype Point Pattern Functions")

# Create test data
create_test_pp3 <- function() {
  set.seed(42)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  coords <- matrix(runif(300), ncol=3)
  marks <- factor(sample(c("A", "B"), 100, replace=TRUE, prob=c(0.6, 0.4)))
  pp3(x=coords[,1], y=coords[,2], z=coords[,3], domain=box, marks=marks)
}

create_single_type_pp3 <- function() {
  set.seed(42)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  coords <- matrix(runif(300), ncol=3)
  marks <- factor(rep("A", 100))
  pp3(x=coords[,1], y=coords[,2], z=coords[,3], domain=box, marks=marks)
}

test_that("is.multitype.pp3 validation works", {
  X_multi <- create_test_pp3()
  X_single <- create_single_type_pp3()
  X_unmarked <- unmark(create_test_pp3())
  
  # Test multitype pattern
  expect_true(is.multitype.pp3(X_multi))
  
  # Test single type pattern
  expect_false(is.multitype.pp3(X_single))
  
  # Test unmarked pattern
  expect_false(is.multitype.pp3(X_unmarked))
  
  # Test invalid input
  expect_error(is.multitype.pp3("not a pp3"))
})

test_that("pp3multitype constructor works", {
  X <- create_test_pp3()
  X_single <- create_single_type_pp3()
  
  # Test with existing marks
  X_multi <- pp3multitype(X)
  expect_s3_class(X_multi, "pp3multitype")
  expect_s3_class(X_multi, "pp3")
  
  # Test with new marks
  new_marks <- factor(sample(c("C", "D"), npoints(X), replace=TRUE))
  X_new <- pp3multitype(unmark(X), marks=new_marks)
  expect_s3_class(X_new, "pp3multitype")
  
  # Test error handling
  expect_error(pp3multitype("not a pp3"))
  expect_error(pp3multitype(X_single)) # Single type should fail
  
  # Test with invalid marks length
  expect_error(pp3multitype(X, marks=factor(c("A", "B")))) # Wrong length
})

test_that("pcfcross3D function interface works", {
  X <- create_test_pp3()
  
  # Test basic function call
  expect_silent({
    g <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=20)
  })
  
  # Check return type
  expect_s3_class(g, "fv")
  expect_true("r" %in% names(g))
  expect_true("trans" %in% names(g) || "iso" %in% names(g))
  
  # Test with different corrections
  g_trans <- pcfcross3D(X, "A", "B", correction="translation", rmax=0.3, nrval=10)
  g_iso <- pcfcross3D(X, "A", "B", correction="isotropic", rmax=0.3, nrval=10)
  
  expect_s3_class(g_trans, "fv")
  expect_s3_class(g_iso, "fv")
  
  # Test error handling
  expect_error(pcfcross3D("not a pp3", "A", "B"))
  expect_error(pcfcross3D(X, "C", "B")) # Invalid mark type
  expect_error(pcfcross3D(X, "A", "C")) # Invalid mark type
  
  # Test with single type pattern
  X_single <- create_single_type_pp3()
  expect_error(pcfcross3D(X_single, "A", "A"))
})

test_that("pcfcross3D produces reasonable values", {
  X <- create_test_pp3()
  
  # Compute pcf
  g <- pcfcross3D(X, "A", "B", rmax=0.4, nrval=20)
  
  # Check basic properties
  expect_true(all(g$r >= 0))
  expect_true(all(is.finite(g$r)))
  
  # Check that theoretical values are 1 (for Poisson)
  if("theo" %in% names(g)) {
    expect_true(all(g$theo == 1))
  }
  
  # Check that computed values are finite
  for(col in names(g)) {
    if(col != "r" && !grepl("theo", col)) {
      expect_true(all(is.finite(g[[col]]) | is.na(g[[col]])))
    }
  }
})

test_that("pcfmulti3D internal function works", {
  X <- create_test_pp3()
  marx <- marks(X)
  
  # Create logical indices
  I <- (marx == "A")
  J <- (marx == "B")
  
  # Test internal function
  expect_silent({
    g <- pcfmulti3D(X, I, J, rmax=0.3, nrval=15, 
                   Iname="Type A", Jname="Type B")
  })
  
  expect_s3_class(g, "fv")
  expect_true("r" %in% names(g))
})

test_that("Edge cases are handled properly", {
  # Test with very small pattern
  set.seed(123)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  coords <- matrix(runif(15), ncol=3) # Only 5 points
  marks <- factor(c("A", "B", "A", "B", "A"))
  X_small <- pp3(x=coords[,1], y=coords[,2], z=coords[,3], domain=box, marks=marks)
  
  # Should work but might produce warnings
  expect_silent({
    g <- pcfcross3D(X_small, "A", "B", rmax=0.2, nrval=10)
  })
  
  # Test with only one point of a type
  marks_rare <- factor(c("A", "A", "A", "A", "B")) # Only one "B"
  X_rare <- pp3(x=coords[,1], y=coords[,2], z=coords[,3], domain=box, marks=marks_rare)
  
  expect_error(pcfcross3D(X_rare, "A", "B")) # Should fail - not enough points
})

test_that("pcfcross3D handles different parameter combinations", {
  X <- create_test_pp3()
  
  # Test with different rmax values
  g1 <- pcfcross3D(X, "A", "B", rmax=0.2, nrval=10)
  g2 <- pcfcross3D(X, "A", "B", rmax=0.4, nrval=10)
  
  expect_s3_class(g1, "fv")
  expect_s3_class(g2, "fv")
  expect_true(max(g1$r) <= 0.2)
  expect_true(max(g2$r) <= 0.4)
  
  # Test with different nrval values
  g3 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=5)
  g4 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=50)
  
  expect_equal(length(g3$r), 5)
  expect_equal(length(g4$r), 50)
  
  # Test with different delta values
  g5 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=20, delta=0.05)
  g6 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=20, delta=0.1)
  
  expect_s3_class(g5, "fv")
  expect_s3_class(g6, "fv")
})

test_that("pcfcross3D handles identical types correctly", {
  X <- create_test_pp3()
  
  # Test with i = j (same type)
  expect_silent({
    g <- pcfcross3D(X, "A", "A", rmax=0.3, nrval=15)
  })
  expect_s3_class(g, "fv")
  
  # Should work but might produce different results than cross-type
  g_AA <- pcfcross3D(X, "A", "A", rmax=0.3, nrval=15)
  g_AB <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=15)
  
  expect_s3_class(g_AA, "fv")
  expect_s3_class(g_AB, "fv")
})

test_that("pcfcross3D handles edge correction methods", {
  X <- create_test_pp3()
  
  # Test translation correction only
  g_trans <- pcfcross3D(X, "A", "B", correction="translation", rmax=0.3, nrval=15)
  expect_s3_class(g_trans, "fv")
  expect_true("trans" %in% names(g_trans))
  
  # Test isotropic correction only  
  g_iso <- pcfcross3D(X, "A", "B", correction="isotropic", rmax=0.3, nrval=15)
  expect_s3_class(g_iso, "fv")
  expect_true("iso" %in% names(g_iso))
  
  # Test both corrections
  g_both <- pcfcross3D(X, "A", "B", correction=c("translation", "isotropic"), rmax=0.3, nrval=15)
  expect_s3_class(g_both, "fv")
  expect_true("trans" %in% names(g_both))
  expect_true("iso" %in% names(g_both))
  
  # Test invalid correction method - should produce a specific error
  expect_error(pcfcross3D(X, "A", "B", correction="invalid"), "unrecognised correction")
})

test_that("pcfcross3D handles different mark types and formats", {
  X <- create_test_pp3()
  
  # Test with character inputs
  g_char <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=10)
  expect_s3_class(g_char, "fv")
  
  # Test with factor inputs (convert to character first to avoid level issues)
  g_fact <- pcfcross3D(X, as.character(factor("A")), as.character(factor("B")), rmax=0.3, nrval=10)
  expect_s3_class(g_fact, "fv")
  
  # Test with numeric inputs (if marks are numeric)
  X_num <- X
  marks(X_num) <- as.numeric(marks(X_num))
  expect_error(pcfcross3D(X_num, 1, 2)) # Should fail - not multitype
})

test_that("pcfcross3D produces consistent results", {
  X <- create_test_pp3()
  
  # Test that multiple calls with same parameters produce same results
  set.seed(123)
  g1 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=20)
  
  set.seed(123)  
  g2 <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=20)
  
  # Should be identical (within floating point precision)
  expect_equal(g1$r, g2$r)
  if("trans" %in% names(g1) && "trans" %in% names(g2)) {
    expect_equal(g1$trans, g2$trans, tolerance = 1e-10)
  }
  if("iso" %in% names(g1) && "iso" %in% names(g2)) {
    expect_equal(g1$iso, g2$iso, tolerance = 1e-10)
  }
})

test_that("pcfcross3D handles domain edge cases", {
  # Create pattern with points near domain boundaries
  set.seed(456)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  
  # Points concentrated near boundaries
  coords <- matrix(c(
    0.1, 0.1, 0.1,    # near min boundary
    0.9, 0.9, 0.9,    # near max boundary  
    0.1, 0.9, 0.5,    # mixed boundaries
    0.9, 0.1, 0.5,    # mixed boundaries
    0.5, 0.5, 0.5     # center
  ), ncol=3, byrow=TRUE)
  
  marks <- factor(c("A", "B", "A", "B", "A"))
  X_boundary <- pp3(x=coords[,1], y=coords[,2], z=coords[,3], domain=box, marks=marks)
  
  # Should work with boundary points
  expect_silent({
    g <- pcfcross3D(X_boundary, "A", "B", rmax=0.5, nrval=10)
  })
  expect_s3_class(g, "fv")
  expect_true(all(is.finite(g$r)))
})

test_that("pcfcross3D handles ratio parameter", {
  X <- create_test_pp3()
  
  # Test with ratio=FALSE (default)
  g_normal <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=15, ratio=FALSE)
  expect_s3_class(g_normal, "fv")
  
  # Test with ratio=TRUE
  g_ratio <- pcfcross3D(X, "A", "B", rmax=0.3, nrval=15, ratio=TRUE)
  expect_s3_class(g_ratio, "fv")
  
  # Ratio version should be a rat object with numerator and denominator attributes
  expect_s3_class(g_ratio, "rat")
  expect_true(!is.null(attr(g_ratio, "numerator")))
  expect_true(!is.null(attr(g_ratio, "denominator")))
})

test_that("pcfcross3D detects clustered species patterns", {
  # Create clustered pattern for species A and B
  set.seed(123)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  
  # Create clusters for species A
  cluster_centers_A <- matrix(c(0.3, 0.3, 0.3, 0.7, 0.7, 0.7), ncol=3, byrow=TRUE)
  coords_A <- rbind(
    cluster_centers_A[1,] + matrix(rnorm(30, sd=0.05), ncol=3),
    cluster_centers_A[2,] + matrix(rnorm(30, sd=0.05), ncol=3)
  )
  
  # Create clusters for species B in similar locations (positive association)
  cluster_centers_B <- matrix(c(0.3, 0.3, 0.3, 0.7, 0.7, 0.7), ncol=3, byrow=TRUE)
  coords_B <- rbind(
    cluster_centers_B[1,] + matrix(rnorm(30, sd=0.05), ncol=3),
    cluster_centers_B[2,] + matrix(rnorm(30, sd=0.05), ncol=3)
  )
  
  # Combine coordinates and marks
  all_coords <- rbind(coords_A, coords_B)
  marks <- factor(c(rep("A", nrow(coords_A)), rep("B", nrow(coords_B))))
  
  X_clustered <- pp3(x=all_coords[,1], y=all_coords[,2], z=all_coords[,3], 
                    domain=box, marks=marks)
  
  # Compute pcf for clustered species
  g_clustered <- pcfcross3D(X_clustered, "A", "B", rmax=0.3, nrval=20)
  
  # For clustered species, pcf should be > 1 at small distances
  if("trans" %in% names(g_clustered)) {
    small_r_indices <- which(g_clustered$r <= 0.1)
    if(length(small_r_indices) > 0) {
      expect_true(any(g_clustered$trans[small_r_indices] > 1.2, na.rm=TRUE))
    }
  }
})

test_that("pcfcross3D detects dispersed species patterns", {
  # Create dispersed pattern for species C and D
  set.seed(456)
  box <- box3(xrange=c(0,1), yrange=c(0,1), zrange=c(0,1))
  
  # Create regular grid for species C
  grid_points <- expand.grid(x=seq(0.1, 0.9, length=3),
                            y=seq(0.1, 0.9, length=3),
                            z=seq(0.1, 0.9, length=3))
  coords_C <- as.matrix(grid_points)
  
  # Create species D points that avoid species C locations (negative association)
  coords_D <- matrix(runif(200*3, 0.1, 0.9), ncol=3)
  # Remove points too close to species C with less aggressive threshold
  dist_to_C <- apply(coords_D, 1, function(pt) min(sqrt(rowSums((coords_C - pt)^2))))
  coords_D <- coords_D[dist_to_C > 0.15, ]
  
  # Ensure we have enough points for both species
  if(nrow(coords_D) < 10) {
    # If filtering removed too many points, generate more
    coords_D <- matrix(runif(300*3, 0.1, 0.9), ncol=3)
    dist_to_C <- apply(coords_D, 1, function(pt) min(sqrt(rowSums((coords_C - pt)^2))))
    coords_D <- coords_D[dist_to_C > 0.1, ]
  }
  
  # Combine coordinates and marks
  all_coords <- rbind(coords_C, coords_D[1:min(20, nrow(coords_D)),])
  marks <- factor(c(rep("C", nrow(coords_C)), rep("D", min(20, nrow(coords_D)))))
  
  X_dispersed <- pp3(x=all_coords[,1], y=all_coords[,2], z=all_coords[,3], 
                    domain=box, marks=marks)
  
  # Compute pcf for dispersed species
  g_dispersed <- pcfcross3D(X_dispersed, "C", "D", rmax=0.3, nrval=20)
  
  # For dispersed species, pcf should be < 1 at small distances
  if("trans" %in% names(g_dispersed)) {
    small_r_indices <- which(g_dispersed$r <= 0.1)
    if(length(small_r_indices) > 0) {
      expect_true(any(g_dispersed$trans[small_r_indices] < 0.8, na.rm=TRUE))
    }
  }
})

# Run all tests
cat("Running 3D multitype point pattern tests...\n")